Analytical solution of a model for complex food webs 
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We investigate numerically and analytically a recently proposed model for food webs [Nature 404, 
180 (2000)] in the limit of large web sizes and sparse interaction matrices. We obtain analytical 
expressions for several quantities with ecological interest, in particular the probability distributions 
for the number of prey and the number of predators. We find that these distributions have fast- 
decaying exponential and Gaussian tails, respectively. We also find that our analytical expressions 
are robust to changes in the details of the model. 
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In ecosystems, species are connected through intricate 
trophic relationships |l], ^| defining complex networks 
H H Hi Hi the so-called food webs. Understanding the 
structure and mechanisms underlying the formation of 
these complex webs is of great importance in ecology 
0i |[ 0- In particular, food web structure provides 
insights into the behavior of ecosystems under pertur- 
bations such as the introduction of new species or the 
extinction of existing species. The nonlinear response 
of the elements composing the network leads to possibly 
catastrophic effects for even small perturbations p| . 

Recently, Williams and Martinez have proposed an ele- 
gant model of food webs — the "niche" model — that just 
with a few ingredients successfully predicts key structural 
properties of the most comprehensive food webs in the 
literature [jlj . Numerical simulations of the niche model 
predict values for many quantities typically used to char- 
acterize empirical food webs that are in agreement with 
measured values for seven webs in a variety of environ- 
ments, including freshwater habitats, marine-freshwater 
interfaces and terrestrial environments. 

We investigate the niche model from a theoretical per- 
spective. We study analytically and numerically the be- 
havior of key quantities for sparse food webs, i.e., webs 
with L <C S 2 , where L is the number of trophic interac- 
tions between species and S is the number of species in 
the web. This is the limit of interest in ecology because 
(i) for most food webs reported in the literature the di- 
rected connectance, defined as C = L/S 2 takes small 
values, and (ii) it corresponds to the limit of large web 
sizes S [|[ ■ We calculate the probability distributions 
of number of prey and of number of predators and find 
that for C <C 1 they depend only on one parameter of 
the model — the average number z of trophic links in the 
network. These distributions give valuable information 
about the structure of the network and enable us 
to calculate other interesting quantities such as the frac- 
tion of "top," "intermediate," and "basal" species, and 
the standard deviation of the "vulnerability" and "gen- 
erality" of the species in the food web H. Our results 
provide compact patterns that describe the structure of 
the food webs generated by the niche model. These pat- 
terns could not have been predicted from the numerical 



simulations reported in Ref. (j| and may be of practical 
and fundamental importance for the study of empirical 
food webs. Moreover, we test our analytical predictions 
with empirical food webs and find agreement. 

Next, we define the model. Consider an ecosystem 
with S species and L trophic interactions between these 
species. These species and interactions define a network 
with S nodes and L directed links. In the model, one 
first randomly assigns S species to "trophic niches" n, 
mapped into the interval [0,1]. For convenience, we will 
assume that the species are ordered according to their 
niche number, i.e., m <ni < ... < ng. 

A species i is characterized by its niche parameter n, 
and by its list of prey. Prey are chosen for all species 
according to the following rule: a species i preys on all 
species j with niche parameters rij inside a segment of 
length Ti centered in a position chosen randomly inside 
the interval [rj/2,nj], with fj = xrii and < x < 1 a 
random variable with probability density function 



p x (x) = b{l-x) 



(6-1) 



(1) 



The values of parameters b and S determine the average 
connectivity z = 2L/S of the food web and the directed 
connectance C = i /S 12 [Jl], [r|. One can also express the 
average number of prey per species as Sr, where the bar 
indicates an average over an ensemble of food webs. It 
then follows that the connectivity is z = 2Sr, the number 
of directed links is L = S 2 r, and the connectance is C = 
r. One can also obtain these expressions in terms of b 
using the equality r — x/2 — 1/ [2(1 + b)]. 

In the niche model, isolated species — that is, species 
with no prey or predators — are eliminated and species 
with the same list of prey and predators — that is, 
trophically- identical species — are "merged" [[hi] . 

First, we address the statistics of the number of prey. 
For large S, the number of prey of a species i is ki = Sri , 
so that the probability distribution p prcy is given directly 
by the distribution of r. Specifically, p pTey (k)=p(r)/S. 

The cumulative probability P(r' > r) = f dr'p(r') is 
the area of the region R of the n — x diagram bounded 
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by lines x = 1, n — 1, and the hyperbole r = nx 

P(r'>r)= I dx ( dn Pn (n)p x (x), (2) 

J r J r fx 

where p n {n) = 1 is the probability density function of n. 
The integration of (||) gives rise to a function involving 
hypergeometric functions |l5[ . To obtain a more "physi- 
cal" solution, one can derivate (H) twice to obtain 



dp(r) 
dr 



Px(r) 



(3) 



In the limit C <C 1, one has b ^> 1 (see [H), so that 
~ be~ bx , and the term in the right-hand side van- 
ishes exponentially, indicating the p(r) and P(r' > r) 
have exponentially decaying tails |]16| . 

To obtain a simpler analytical solution for p(r) than 
given by the hypergeometric functions, we approximate 
p x in the entire arrange by an exponential. We expect 
the results to be the same for x <C 1 (U because p x 
takes non- vanishing values only for small x. Under this 
approximation, the integration of (j^) yields 



p(r) =bEi(br), 



(4) 



where Ei(x) — f dt t~ x exp(— t) is the exponential- 
integral function |D| . The probability distribution 
Pprcy(k) is obtained from (^) making the substitutions 
r = k/S and b = S/z, yielding 

PpIey (k) = {l/z)E 1 {k/z). (5) 

We compare in Figs. 0(a) and (b) the predictions of (^) 
with numerical simulations. We find close agreement be- 
tween our analytical expression and the numerical re- 
sults. In particular, they show an exponential decay for 
large k. The deviations observed for small values of k 
are due to the fact that kj — Srj is an average value 
implying that it is a good approximation only when the 
fluctuations of kj are small, which is no longer true for 
small k. 

Next, we address the statistics of the number of preda- 
tors. Note that for f <C 1 (li), the predators of species 
i have, to first approximation, niche values rij > rij and 
that the segment rj is placed with equal probability in the 
interval [0, rij]. Therefore, the probability for a species j 
to prey on i is rj/nj = XjUj/rij — Xj, implying that the 
average probability for the species with rij > rij to prey 
on species i is x. 

If we assume that S 3> 1, the number of predators of 
i is the result of S — i independent "coin-throws" with 
probability x of being a predator and probability 1—3? 
of not being a predator, implying that the probability of 
species i having m predators is given by the binomial dis- 
tribution. It then follows that the distribution of number 
of predators for a general species is the average over the 
different binomials 
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FIG. 1: (a) Linear and (b) log-linear plots of the distribution 
of number of prey for 1000 simulations of food webs with S = 
1000. We show results for z = 10, 20 and the corresponding 
theoretical predictions. As expected, we find an exponential 
decay of the distributions, (c) Linear and (d) log-linear plots 
of the distribution of number of predators for the same food 
webs as in (a) and (b). As predicted, we find a regime where 
the distribution is uniform followed by a Gaussian decay. We 
test our analytical predictions with empirical data |l| for (e) 
Pprey(k) and (f) Ppred(w-) for Bridge Brook (solid line) and St. 
Martin (broken line). 



In the limit of interest, S > 1, S< 1, and Sx = z, one 
can approximate the binomial distribution by a Poisson 
and the sum by an integral 



Pprcd(m) = - 

z 



dt 



f 



~7(m+l,*), (7) 
z 



where 7 is the "incomplete gamma function" flTa] , |l7fl . For 
m < z/2, the function 7 is approximately constant, while 
it decays with a Gaussian tail for m ~ z. In Fig. |l|(c) 
and (d), we compare the predictions of (0) with numerical 
solutions and find good agreement. 

In Figs. |l|e,f, we compare our analytical predictions, 
Eqs. (H)-(0), with data from recent food webs: Bridge 
Brook (S = 25, z = 8.6) and St. Martin Island (S = 
42, z = 9.8). We find that the distributions of number of 
prey is well approximated by the data and that the distri- 
butions of number of predators are "noisy" but still show 
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the expected cutoff for m ~ z and is approximately con- 
stant for m < z as predicted by Eq. (0). This agreement 
is remarkable since the webs analyzed are quite small so 
one might not expect the theoretical expressions to hold. 

Next, we evaluate the fraction of top T, intermedi- 
ate / and basal B species. As the names indicate, top 
species have no predators and basal species have no prey, 
while intermediate species are those with both prey and 
predators. The fraction of intermediate species is just 
I = 1 — (T + B). The fraction T of top species is, by 
definition, 

1 — exp(— z) 



Ppred(O) 



(8) 



Note that a similar result is obtained if one calculates the 
sum (||) for m — 0. Since typically 5 < z < 20, Eq. (||) 
can be approximated simply as T = 1/z. 

To calculate the fraction B of basal species, we note 
that a species has no prey only if its range r falls in a 
region with no species |L8|. In the limit of large S, the 
probability density for finding an empty interval of length 
<5 is Se~ ss , as predicted by the canonical distribution p9[ . 
Thus, the probability of finding a species-free segment of 
length larger than r is e~ Sr , which gives the probability 
for a species of range r not to prey on other species. 
Using (ij), it follows that the average probability is: 
-Sr ^ Wl + z) 



B = 



dr i 



Mr) 



(9) 



In the model Q], isolated species are eliminated, so 
they are not counted towards top or basal species. To 
correct the estimates (||)-(||) for this effect, we remove 
the isolated species. We estimate the number of isolated 
species to first order by assuming that having no prey is 
statistically independent of having no predators, imply- 
ing that the fraction of isolated species is just the product 
of the fractions of top and basal species. This assump- 
tion does not take into account the possibility that a 
species with no prey is likely to have a low niche value 
n and hence it has a high probability to have predators. 
Nonetheless, this simple approximation provides an up- 
per bound for the number of isolated species which leads 
to a lower bound on T and B, 

. T-TB . B-TB 

T = , B 1 = (10) 

1 — TB l-TB v ' 

In Fig. [2J we compare our analytical predictions for the 
fraction of top and basal species with numerical simula- 
tions of the model. As expected, Eqs. (^)-(|To|) provide 
bounds for the numerical results. 

Finally, we calculate the standard deviations of the 
vulnerability and generality of the species in food webs 
generated according to the model. The vulnerability of a 
prey is defined as its number m of predators, and the gen- 
erality of a predator as its number k of prey. Following 
Rcf. [Q, we define the normalized standard deviations of 

the vulnerability as a v — m 2 /m 2 — 1 and of the generality 

2 

as <Jq = k 2 /k — 1. By definition, one has m = k = z/2 

for both cases. 
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FIG. 2: Fraction of top and basal species as a function of 
the average connectivity z. The shaded region corresponds to 
the interval of z typically observed in empirical food webs, (a) 
Comparison of the results of 100 simulations of food webs with 
S — 1000 — for which isolated species where not removed — 
with our theoretical predictions Eqs. Note the good 

agreement between the analytical expressions and the numer- 
ical results, (b) Comparison of the results of 100 simula- 
tions of food webs with S = 1000 — for which isolated species 
were removed — with our theoretical predictions Eqs. di|)-(p^|). 
Note that the theoretical predictions provide narrow bounds 
for the numerical results. 



To evaluate ay, we first calculate m 2 . Equation (0) 
yields ni 2 — z 2 /3 + z/2, so that 



(11) 



We next calculate ctg, for any value of C, by direct 
evaluation of k 2 . If 5 > 1, the number of prey of a 
species having a range r is k = Sr, and we find that 

F/fc 2 = = 8(6 + l)/[3(6 + 2)], implying that 



S 



1 



31 + 2C 



- 1. 



(12) 



For C C 1, (Tg becomes a constant with value y5/3, 
a result that can also be obtained from (g). We show 
in Fi g. j3| the results for our analytical expressions (|ll| ) 
and (]l2|) and compare them with results from numerical 
simulations of the niche model. 

We have also studied the robustness of our predictions 
to changes in the particular formulation of the details 
of the model. The nature of approximations used in 
the derivations of the expressions for the distributions 
of the number of prey and predators, Eqs. (5)- (7), allow 
us to conclude that: (i) the distribution of the number of 
predators does not depend on the specific form of p(x) . 
The only requirement is that the connectance C = x/2 
tends to zero under some limit, so that z = SC remains 
finite when S tends to infinity; and (ii) the distribution 
of the number of prey depends on the functional form of 
p(x), but Eq. (7) will still be is obtained for all p(x) de- 
caying exponentially as C tends to zero. Thus, it appears 
that our findings are robust under quite general condi- 
tions, a result that is not possible to obtain without an 
analytic treatment of the problem. 
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FIG. 3: Normalized standard deviations of generality and vul- 
nerability as a function of the average connectivity z. The 
shaded region corresponds to the interval of z typically ob- 
served in empirical food webs, (a) Comparison of the re- 
sults of 100 simulations of food webs with S = 1000 — for 
which isolated species where not removed — with our theoret- 
ical predictions Eqs. (|ll"|)-(|l2]). (b) Comparison of the results 
of 100 simulations of food webs with S = 1000 — for which 
isolated species where removed — with our theoretical predic- 
tions Eqs. Note that as for Fig. H, removing iso- 
lated species leads to slightly less good agreement with the 
simulation results for cry . However, the removal of isolated 
species does not appear to be a factor in the deviations found 
for oq. The reason why o~c underestimates the simulation 
results at small z values relates to the fact that kj = Srj is 
a good approximation only when the fluctuations of kj are 
small, which is no longer true for small k. 



Our results are also of interest for a number of other 
reasons. First, we demonstrate for the first time that the 
distributions of number of prey and number of predators 
have different functional forms. Second, we show that 
both distributions have characteristic scales, i.e., both 
distributions have well defined means and standard devi- 
ations as S increases to infinity. Third, we find that the 
functional forms of the distributions of number of prey 
and number of predators depend only on the average con- 
nectivity z, and agree with empirical data. This result is 
rather surprising in face of the complexity of the empir- 
ical and model food webs. Finally, we show that other 
quantities of biological interest also depend exclusively z. 
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